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ABSTRACT 

The  general  3-D  inversion  formula  developed  by  Bleistein  et  al  is  specialized 
to  a  depth-dependent  medium.  An  efficient  ray  tracer  for  this  model  is  developed 
to  calculate  the  necessary  constituents  of  the  ray  theory.  The  traveltime  and 
amplitude  tables  for  the  inversion  are  computed.  The  motivating  application  for 
this  project  is  detection  of  small  objects-lOcm-in  shallow  water-  10-20m. 


INTRODUCTION 

This  paper  describes  the  background  mathematical  analysis  for  specialization  of 
the  Born/Kirchhoff  inversion  formalism  [Cohen,  Hagin  and  Bleistein,  1986;  Bleistein, 
1986],  to  a  three  dimensional  (3-D)  depth  dependent  background  medium.  Some 
numerical  results  that  test  the  code  I  developed  to  implement  this  method  are  also 
presented. 

The  underlying  motivating  project  has  as  its  objective  the  detection  of  small  ob¬ 
jects  (10-20cm)  in  shallow  water  (10-20m).  This  objective  differs  from  the  seismic 
inverse  problem  in  two  ways:  first,  the  depth,  time  and  frequencies  differ  from  the 
seismic  problem  by  approximately  a  factor  of  100 — not  an  important  difference;  sec¬ 
ond,  the  lateral  extent  of  the  scatterer  is  finite  and  small  compared  to  the  range  from 
the  source/receiver  array  or  survey  extent,  whereas  in  the  seismic  problem,  the  lateral 
extent  is  often  comparable  to  these  other  length  scales.  We  have  not,  at  this  time, 
exploited  this  latter  difference,  so  that  the  resulting  code  is  applicable  to  both  types 
of  problems  at  both  scales. 

The  survey  we  have  in  mind  is  traditional:  a  boat  carries  a  towed  array  along 
parallel  lines  on  the  upper  surface  and  periodically  sets  off  an  acoustic  source,  collect¬ 
ing  the  resulting  backscattered  data  on  the  array.  Data  is  then  re-sorted  in  common 
(constant)  offset  data  sets  and  processed  by  the  program  developed  here  to  produce 
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a  reflector  map.  In  addition,  the  output  provides  an  estimate  of  the  angularly  depen¬ 
dent  reflection  coefficient  at  specular  at  the  sample  points  on  the  scatterer,  as  well  as 
aji  estimate  of  that  incident  specular  angle  with  respect  to  the  normal  to  the  reflector. 
By  processing  the  data  for  a  suite  of  offsets,  data  for  amplitude  versus  offset  (AVO) 
or  amplitude  versus  angle  (AVA)  analysis  is  generated  for  each  point  on  the  reflector. 

To  calculate  the  amplitude  weights  for  a  3-D  common  offset  v(z)  inversion  pro¬ 
gram,  I  use  the  analytic  formulas  for  3-D  ray  data  in  a  depth  dependent  propagation 
medium.  These  ray  data  are  interpolated  from  ray  coordinates  to  cartesian  coordi¬ 
nates  after  ray  tracing.  From  these  ray  data,  I  compute  the  Jacobian  of  the  transfor¬ 
mation  from  Cartesian  {x,  y,  z)  coordinates  to  ray  coordinates  (<7,  a,  (5)  and  I  compute 
the  Beylkin  determinant  h.  Both  of  these  determinants  are  essential  to  the  amplitude 
computation. 

The  numerical  tests  confirm  that  the  computer  code  matches  analytical  values 
of  the  travel  time  and  amplitude  along  the  rays  extremely  well.  Furthermore,  for 
a  model  sphere  in  constant  background,  we  find  that  the  reflector  map  successfully 
produces  an  image  of  the  sphere. 


INVERSION  FORMULA 


In  this  section,  the  inversion  formulas  of  Bleistein  [1986]  are  introduced.  These 
formulas  provide  a  tool  for  obtaining  correct  locatios  of  interfaces  as  well  as  a  model- 
consistent  specular  reflection  coefficient  and  incidence  angle.  The  inversion  formulas 
have  the  form  of  aperture-limited  Fourier-like  integrals.  The  integrand  of  these  inte¬ 
grals  contains  a  determinant  that  characterizes  the  viability  of  inverting  a  particular 
data  set.  This  determinant  is  part  of  a  Jacobian  that  depends  on  both  the  back¬ 
ground  propagation  parameters  and  on  the  source-receiver  configuration.  As  a  result, 
the  problem  of  extending  the  inversion  formula  to  new  recording  geometries  is  reduced 
to  a  problem  of  computing  the  value  of  the  form  of  this  determinant  that  corresponds 
to  a  specific  geometry. 

The  full  3-D  Kirchhoff  inversion  formula  Bleistein  [1986]  is 


^a(y,O|V,0(y,OI 


j  iu>  du!  e  Ug {Xg ,Xs,u). 


(1) 


In  this  equation,  /?(y)  is  the  reflectivity  function  for  the  imaging  section;  y  =  (j/i,  ^2,  Vs) 
is  the  3-D  position  vector  of  the  output  or  imaged  point  and  ^  =  (^1,^2, 0)  are  the 
surface  coordinates  of  the  ray  to  y  from  either  the  source  or  receiver.  For  a  common 
offset  inversion,  we  denote  by  As  and  {Ag  and  Tg)  the  amplitude  and  phase  of  the 
WKBJ  Green’s  functions  at  y  with  initial  point  at  the  source.  Thus,  in  the  above 
equation 

a(y ,  0  =  ^sAg,  0(y,  0  =  'T,  -F  Tg.  (2) 

with  the  latter  being  the  total  traveltime  from  source  to  y  to  receiver.  Furthermore, 
Us{Xg,Xs,u!)  is  the  data  in  u>  (frequency)  domain;  and  Xg  are  the  source  and 
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the  receiver  position  coordinates;  h  is  the  determinant  characterize  the  viability  of 
inverting  a  particular  data  set  with  a  specific  geometry  at  the  given  output  point.  The 
expression  of  the  spatial  weighting  in  terms  of  this  one  determinant  for  any  source- 
receiver  configuration  and  background  propagation  speed  is  a  major  contribution  of 
Beylkin’s  approach  to  high-frequency  inversion  and  is  being  referred  as  the  ’’Beylkin 
determinant” .  The  details  of  h  will  be  described  later. 

The  reflectivity  function  provides  a  reflectivity  map  through  a  family  of  bandlim- 
ited  delta  functions  that  peak  on  the  reflectors  in  the  medium.  The  peak  amplitude 
is  proportional  to  the  angularly  dependent  reflection  coefficient  at  a  specular  angle 
Os,  where  29s  is  the  angle  between  the  ray  directions  from  the  source  and  receiver  to 
the  output  point  y: 

Ppeak--  R{y,9s)—^  I  F{uj)d0J.  (3) 

TTCfy)  7-00 

That  is,  the  peak  value  of  (3  for  y  on  S'  is  the  geometrical  optics  reflection  coefficient 
multiplied  by  2cos0s/c(y)  and  multiplied  by  l/27r  times  the  area  under  the  filter  in 
the  w-domain. 


|V,./.(y,OI 


For  y  not  on  the  reflector  surface,  the  angle  9  is  defined  by  the  equation, 

V.r(y.x,).V,r(y,x,)  =  ^.  (4) 

Then,  by  computing  V0(y,^)  •  V(/)(y,^),  we  are  able  to  show  that 

.  ,  . .  2  cos  9  ,  , 

(5) 

In  particular,  for  y  on  the  reflector  surface,  this  provides  a  means  for  estimating 
cos  ^3,  through  the  introduction  of  another  inversion  operator,  differing  from  (3  by 
one  power  of  |Vj^0|: 

=  i^/'''^a(y.O|v;^L)P  (6) 

The  addition  of  the  extra  divisor  of  \Xy(t>\  introduces  this  divisor  in  the  asymptotic 
amplitudes  of  the  result.  In  particular; 


J  iuj  du)  e 


Plpeakiy)  ~  R{y,9s)-^  J  F{uj)  du, 


The  fact  that  these  two  operators  differ  by  a  factor  of  2cos9/c{y)  allows  us  to  esti¬ 
mate  cos  9s  from  the  ratio  of  the  outputs  without  ever  having  determined  the  specular 
source- receiver  pair  that  produced  the  distinguished  value  of  9s.  This,  in  turn  allows 
us  to  determine  R{y,9s)  from  either  output.  With  knowledge  of  9s  and  the  back¬ 
ground  wavespeed,  c(y),  it  is  conceivable,  within  the  limits  of  the  accuracy  of  the 
data,  that  we  should  be  able  to  estimate  the  jump  in  the  propagation  speed  across 
the  reflector  using  the  formula  for  the  geometrical  optics  reflection  coefficient.  More 
generally,  in  a  variable  density  medium,  we  determine  the  impedance  jump  from  this 
analysis.  By  processing  data  for  multiple  offsets,  we  generate  data  for  amplitude 
versus  offset  (AVO)  analysis  or  amplitude  versus  angle  (AVA)  analysis. 
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Fig.  1.  Diagram  of  coordinate  system  showing  an  initial  ray  direction. 

RAY  TRACING 

Here,  I  describe  the  ray  tracing  in  a  depth-dependent  medium  in  3-D.  For  a  more 
complete  discussion  of  the  underlying  ray  theory  for  determining  solutions  of  the  wave 
equation  in  the  form  Aexpfio^r],  refer  to  Bleistein  [1984]. 


The  ray  equations 

The  general  form  of  the  ray  equations 


dxi 
da 
dr 

Ta=P 


=  Pi, 
2 


dpi  _  dp 
da  dxi  ’ 


i  =  1, 2, 3, 


+P2^  +Ps^ 


(8) 


provide  the  basis  for  ray  theoretic  modeling. 

The  slowness  vector  points  in  the  direction  normal  to  the  surfaces  of 

constant  r.  Surfaces  of  constant  r  are  called  “wavefronts” ,  and  the  p  vector  points  in 
the  direction  tangent  to  the  “raypaths” .  These,  in  turn,  are  the  spatial  trajectories 
of  the  solution  to  (8). 

We  seek  solutions  of  (8)  for  rays  emanating  from  a  single  point,  say  aro,  yo,  zq  in  an 
arbitrary  downward  direction.  Those  directions  are  determined  by  the  initial  values 
of  a  and  (3  representing  the  azimuthal  and  polar  angles,  respectively,  of  the  spherical 
coordinate  system.  See  Fig.  1.  Then,  the  initial  ray  direction  po  is  given  by 


Po  =  po(cosasin/?,sinasin/?,cos/?). 


(9) 
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Specialization  to  a  Depth-dependent  medium 

For  a  medium  that  has  wavespeed  variability  only  in  the  z  direction,  pi  and  p2 
are  constants  on  each  ray  path;  denote  them  by  pio  and  p20-  Note  that  pio  and  p20 
are  their  respective  initial  values. The  ray  equations  become 


dxi 

da 


=  Pio, 


dpi 

da 


=  0, 


dpz  _  d£_ 

da  ^  dx^' 


dr 

da 


From  the  eikonal  equation,  we  know  the  value  of  pz 


P3 


\Jp^{xz)  -Pio  -P20 
\Jp^{xz)  -p‘^{xzQ)s\n^  (5. 


(10) 


(11) 


The  third  equation  in  (10)  relates  a  and  xz,  so  the  ray  equations  may  be  rewritten  in 
terms  of  X3: 


dxi  _  dxz  _  P20 

dxz  pz  ’  dxz  pz  ’ 

clx^  dx^ 

dpz  _  p  dp 
dxz  Pz  dxz ' 
dr  p^  da  1 

dxz  Pz '  dxz  Pz ' 


(12) 


We  solve  for  xi,  xz,  t,  and  a  as  functions  of  xz,  o,  and  /?.  Let  the  initial  depth 
be  X30  =  zo  and  the  final  depth  be  xz  =  z.  The  solution  of  (13),  then,  is 


X  —  Xq 


y-vo 


cos  a  sin  /?  dz' 

c{zq)  Jzo  ^p’^[z')  —  p‘^(zo)  sin^  /? 

sin  a  sin/?  dz' 

c{zq)  Jzo  ^p^(^z')  —  p‘^{zq)  sin^  /? 


r-To 


_ ^ _ 

{z')^p‘^{z')  -p2(zo)sin^^’ 


(13) 
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dz' 

cr  —  (To  =  /  ,  „  ••• 

\Jp^{z')  -p^{zq)s\v?  (3 

Recall  that  a  and  /?  are  the  initial  angles  of  the  ray.  When  pz  =  0,  the  ray  is 
horizontally  propagating.  This  is  called  a  “turning  point” .  The  equations  above  all 
have  integrable  singularities  at  the  turning  point,  requiring  special  care  in  numerical 
computation,  to  be  discussed  below. 


The  Jacobian  J  and  Ray  Amplitude 

Here,  I  discuss  the  solution  of  the  transport  equation 

2Vt  •  VA  +  AVV  =  0,  (14) 

for  the  amplitude  A.  This  equation  can  also  be  written  as  an  ordinary  differential 
equation  in  the  ray  parameter  a 

=  0,  (15) 

do 


with  solution  (Bleistein  1984,  8.3.12) 

A  = 

47r^c(0J(<7,a:,/?) 


(16) 


In  these  equations,  J  is  the  Jacobian  of  the  transformation  from  x  to  (cr,  a,  /?)  via  the 
solution  of  the  ray  equations. 


J  = 


d{o,  a,  P) 


(17) 


We  calculate  the  ray  theoretic  amplitude  A  by  computing  J. 

Let  us  consider  the  solution  of  the  ray  equations  in  terms  of  cr,  rather  than  2:. 
The  solution  is  a  family  of  rays,  distinguished  from  one  another  by  the  choice  of  the 
parameters  a  and  p.  Along  each  ray,  the  values  of  r,  cr,  x,  y,  z,  p,-  are  known.  pi  and 
P2  are  constants  on  each  ray;  that  is,  independent  of  cr: 


Pi 


— — r  COS  OL  sin  /?, 
c{Zq) 


P2  =  — r  sin  o;  sin/?. 

c{zo) 


(18) 
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With  Pi  and  p2  independent  of  cr,  the  equations  for  x  and  y  are 


a:  =  ^1  H — ; — r  cos  a  sin  /?, 
c{zo) 


y  -  ^2  +  — r  sin  a  sin  /?. 

c{zo) 


(19) 


and  2  is  given  in  terms  of  a  by 

z  =  j  Pzda.  (20) 

I  do  not  use  equation  (20)  since  we  do  the  calculation  on  a  uniform  2  grid,  instead, 
we  calculate  a  for  each  z  by  using  the  third  equation  in  (10).  From  (20)  and  (10),  we 
can  derive  the  nine  terms  of  the  determinant  J  in  (17)  as 


dx 

da 

sin  /?  cos  a 
c{zo)  ’ 
sin/?  sin  q; 

da 

dx 

c{zo) 

dz 

Ta=^^' 
asm /3  sin  a 

da 

c(^o)  ’ 

^  _ 

(7sin  /?  cos  a 

da 

c(zo)  ’ 

dz 

w-  =  0, 

da 

dx 

acos  (3  cos  a 

c{zo) 

dy 

acos  /?  sin  Q 

dd~ 

c(zo)  ’ 

dz 

f  —  sin  /?  , 

/  ,  .da. 

J  c{zo) 

(21) 


With  these  values,  we  can  calculate  J  in  (16)  and  then  calculate  A  by  using  (17). 
Thus,  we  obtain  A  on  a  uniform  grid  in  2,  but  calculate  it  as  if  the  independent 
parameter  along  the  ray  were  o. 


The  Belkyin  determinant  h 

The  effects  of  source-receiver  geometry  are  completely  described  by  the  Beylkin 
determinant 
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Ps  +  P5 


h{y^O  =  det 


^(P.  +  PJ 

56 


(22) 


5(P.  +  Pq) 

-  56  - 

in  terms  of  the  slowness  vectors  and  and  their  derivatives  with  respect  to  the 
surface  coordinates,  (6)6)-  These  vector  quantities  are  evaluated  at  the  output  point 


y- 

The  determinant  must  be  finite  and  nonzero  for  the  identification  of  the  cascaded 
model  and  inversion  integral  as  an  approximate  Fourier  transform.  Thus,  we  could 
use  the  value  of  this  matrix  to  characterize  source-receiver  configurations  as  providing 
invertible  data  by  this  formalism  at  an  output  point,  y.  In  particular,  we  require  that 
this  determinant  be  finite  and  nonzero  for  some  range  of  ^  values  at  any  y  where  the 
high  frequency  inversion  is  to  be  computed. 

The  general  form  of  the  Beylkin  determinant  can  be  written  for  all  source  receiver 
geometries  as 


Hy,  0  =  {Ps  +  Pg) '  +  '^9)  ^  +  '^9) 

=  {Ps  +  Pg)  -VgXWg-i- 

{Ps  +  Pg)  -  Vs  XWs-\-  (23) 

{Ps  +  Pg)  -VgXWs-i- 
{Ps  +  Pg)  •  Vs  X  Wg. 


where. 


Vs  = 


dps 

56’ 


dPg 


"""  "  56  ’ 


dp. 

w.  = 


=  ^ 
56’  "“56' 


(24) 


The  first  two  terms  can  be  recognized  as  being  the  Beylkin  determinants  for  the  com¬ 
mon  shot  and  common  receiver  geometries.  We  may  write  the  Beylkin  determinant 
as 


h{y,  0=2  cos^  e  [hs{y,  0  +  0]  +  (25) 

{Ps  +  Pg)  •  XWs-\-VsX  Wg]  , 

where  hg{y,^)=Pg  ■  Vg  x  Wg  is  the  determinant  for  the  common  receiver  case  and 
hs{y,0—Ps  •  Vs  X  Ws  is  the  determinant  for  the  common  source  case.  Unfortunately, 
the  last  term  in  (26)  contains  the  cross  products  Vg  x  Wg  and  Vg  x  Wg,  which  are  not 
easily  simplified. 
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Another  way  of  writing  the  Beylkin  determinant  is 

h(y,0  =  2cos^0lh,(y,^)  -hhgiy,^)] 
+ha(p,0  +  hb(p,0 
0  d"  hdivj 0> 

where  the  ha,  hb,  he  and  hd  are  defined  by 

ha{y,0  =Ps-VsX  Wg,  hb{y,0  =Ps-VgX  Ws, 
hc{y, 0  =  Ps  •  X  hd{y,  =pg- Vs  X  Wg. 


(26) 


(27) 


Unfortunately,  we  cannot  write  the  general  common-offset  Beylkin  determinant  in  a 
form  that  has  a  common  multiplier  of  cos^  9. 

We  already  know  the  vectors  and  Pg  at  the  output  point  from  the  solution  of 
the  ray  equations.  Now  we  will  compute  their  partial  derivatives  with  respect  to 
and  ^2  as  solutions  of  another  set  of  ray  equations. 

By  taking  derivatives  of  (19)  with  respect  to  and  ^2>  respectively,  and  using 
(20),  we  obtain  differential  equations  for  dpifd^,  and  dp2ld^i.  These  equations  are 


dpi 


dp2 

di2 


dp2  do  . 

%  = 

dpi  do 

'XT  - 

56  56 

do  . 


(28) 


The  only  quantities  that  we  do  not  know  in  the  above  equations  are  dofd^i 
and  dofd^2-  It  is  difficult  to  implement  the  derivatives  with  respect  to  6  and  6- 
However,  using  the  symmetry  of  the  v{z)  medium,  we  find  that  d/d^i  =  —dfdx  and 
d/d^2  =  —dfdy  as  shown  schematically  in  Fig.  2.  That  is,  when  a  ray  is  shifted 
horizontally,  the  quantity  o  on  that  ray  does  not  change.  Also,  one  can  see  from 
Fig.  2,  dofd^  =  \oi  —  02)/d^  =  — (cr2  —  oi)ldx.  Thus,  dofd/^i  and  5(t/9/6  are 
determined  from  results  already  calculated. 

From  the  eikonal  equation-the  last  equation  in  (8),  the  derivatives  of  with 
respect  to  6  and  6  can  also  be  found: 


5p3 

56 


5pi  5p2w 
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Fig.  2.  Ray  symmetry  used  in  calculating  h 


dp3 


dpi  dp2 


The  derivatives  for  source  and  receiver  are  calculated  seperately  using  (28)  and  (29). 
There  is  a  certain  symmetry  here  that  can  be  exploited  since  the  difference  between 
the  source  and  the  receiver  is  just  a  shift  in  x.  Once  we  have  h,  together  with  the 
product  of  ray  theoretic  amplitudes,  a  in  1  and  2,  the  amplitude  computation  is 
complete. 


NUMERICAL  IMPLEMENTATION 

Ray  data  are  generated  along  each  raypath  with  a  unique  coordinate  reference 
(a,  /?,  a).  Transformation  of  ray  data  from  this  ray  system  to  a  uniform  grid  is  achieved 
by  linear  interpolation  as  shown  in  Fig.  3.  Rays  intersect  with  a  z-plane.  Four 
adjacent  intersection  points  that  surround  a  grid  point  are  found.  Linear  interpolation 
is  used  to  determine  the  ray  data  at  the  grid  point  from  the  four  intersection  points. 

This  scheme  is  not  accurate  in  the  vicinity  where  rays  are  near  turning  because 
the  distance  between  the  intersection  points  will  be  too  large.  To  deal  with  rays  prop¬ 
agating  nearly  horizontally,  that  is,  when  pz  is  small,  vertical  interpolation  is  used. 
Vertical  interpolation  differs  from  the  above  scheme  in  that  it  use  a  vertical  plane  in¬ 
stead  of  a  horizontal  2:-plane.  The  nearly  horizontally  propagating  rays  will  intersect 
the  vertical  plane  with  a  closer  spacing  leading  to  a  more  accurate  interpolation  than 
if  we  used  a  horizontal  interpolation  grid. 

Now  we  have  all  the  ray  data  we  need  on  the  uniform  {x,y,z)  grid.  Use  these 
ray  data  to  calculate  the  Jacobian  J,  ray  amplitude  a  and  Beylkin  determinant  h. 
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surface 


Combining  these  values  gives  the  amplitude  weight  value  at  each  grid  point  for  that 
particular  source  point.  This  amplitude  table  along  with  the  traveltime  table  will  be 
used  in  the  inversion,  it  is  just  a  shift  and  interpolation  for  other  source  or  receiver 
points  as  the  ray  data  are  invariant  for  lateral  displacement. 


EXAMPLE:  CONSTANT  WAVESPEED 


In  this  section,  we  present  the  specialization  of  the  above  results  for  a  constant 
wavespeed  medium  where  all  calculations  can  be  carried  out  explicitly,  in  this  case, 
equation  (17)  reduces  to 

J  =  —  sin  (3,  (30) 

c 

and  equation  (16)  becomes, 

A  =  — . 

47rr 

For  the  special  case  of  constant- wavespeed,  with  a  generally  non-fiat  surface,  some 
simplification  of  h  is  possible.  Application  of  Gaussian  elimination  produces  the 
following  simplified  forms  for  ha,  hi,,  he,  and  h^  in  (27) 


ha{y,0 

hb{y,0 

hc{y,0 


c^s  '  ■  56  56  ’ 

-1  ..  dpg  dXs 

56  56 

-1  ,  dps  dxg 

ch'g^^  56  ^  56’ 


(32) 
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.  /  -1  - 

For  the  special  case  of  constant  wavespeed,  with  a  flat  recording  surface,  hg,  hg,  ha, 
hb,  he,  and  hd  become 


hn  - 


Vs 


hg  = 


2/3 


u  - 
h„  = 


a  —  Q  o  ’ 
(yrjTg 


_  t/3Cos2^ 
hb  = 


(33) 


_  ^3  cos  20  u  - 

-  or,  1  - 


c^r'^Tg 


c?rgrl 


Substituting  these  results  into  equation  (26),  the  full  expression  for  the  Beylkin  de¬ 
terminant  may  be  written  as 


^(2/,0  =  2cos2  0  ^ 


(^’^  + 


'  s*  g 


(34) 


The  corresponding  formulas  for  /?(?/)  and  (5i{y)  for  3D  common-offset,  constant- 
wavespeed,  and  with  a  flat  recording  surface  become 


{rs  +  rg){rj  +  rj) 


j.2-2 

>9 


COS  9 


j  iu  du  e  '^^'^‘'^''^^^"'us{xg,Xa,u!), 


(35) 


and 


{ts  +  rg){rl  +  r^) 


'  s'  g 


J  iu  duj  e  "^^''‘'^''<'^^‘'us{xg,Xs,u). 


(36) 


Equation  (36)  exactly  matches  equation  (30)  in  Sullivan  and  Cohen  [1987],  which  was 
derived  using  a  different  approach. 

Below  is  a  list  of  comparison  of  the  numerical  results  and  the  analytic  results  of  a 
constant  background.  The  grids  are  40  x  40  x  40  and  the  ray  shooting  grids  in  a  and 
P  are  15  x  15,  which  is  rather  coarse.  The  numerical  results  have  two  or  three  digit 
accuracy: 
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Depth 

Analytic 

Numerical 

Error 

4A^ 

0.059402 

0.062488 

5.19% 

8A2 

0.287655 

0.291268 

1.26% 

12Az 

0.483947 

0.490268 

1.31% 

16A^ 

0.638489 

0.644558 

0.95% 

20A2 

0.754078 

0.759777 

0.76% 

24A2 

0.838724 

0.843727 

0.60% 

28A2 

0.900615 

0.905202 

0.51% 

32A^ 

0.946287 

0.950361 

0.43% 

36A^ 

0.980478 

0.983688 

0.33% 

40Az 

1.017178 

1.020043 

0.28% 

These  are  the  amplitude  weights  on  one  trace  (with  a  4-sample  skip)  at  offset 
10 A2:.  Notice  that  the  accuracy  increases  with  depth,  this  is  because  the  polar  angle 
P  associated  with  the  rays  increases  with  depth  and  thus  has  a  closer  ray  spacing  for 
interpolation.  The  computation  of  travel  time  was  about  10  times  more  accurate. 

INVERSION  EXAMPLES 

The  first  example  provides  a  simple  test  of  positioning  and  amplitude  accuracy. 
The  model  consists  of  constant  velocity  layers,  as  shown  in  Figure  4.  Layer  velocities 
are  2,  3,  6  and  lOkm/s  from  the  top  to  the  bottom.  Figure  5  shows  the  model  data 
with  oflFset=3km.  The  inversion  produces  two  outputs:  one  for  and  one  for  Pi. 
The  amplitudes  in  the  two  outputs  differ  by  the  factor  cos  9,  where  9  is  the  specular 
angle.  For  the  first  interface,  the  real  values:  R  is  0.200,  Rcos0  is  0.172  and  9  is  30.90 
degrees.  The  inversion  results:  R  is  0.203,  Rcos^  is  0.174  and  9  is  30.84  degrees.  The 
percentage  errors  are  1.5,  1.1  and  0.2  percent  respectively.  For  the  second  interface, 
the  exact  values:  R  is  0.333,  Rcos^  is  0.307  and  9  is  22.76  degrees.  The  inversion 
results:  R  is  0.341,  Rcos0  is  0.315  and  9  is  22.71  degrees.  The  percentage  errors  are 
2.4,  2.6  and  0.2  percent  respectively.  For  the  third  interface,  the  exact  values:  R  is 
0.250,  Rcos^  is  0.235  and  9  is  19.94  degrees.  The  inversion  results:  R  is  0.256,  Rcos^ 
is  0.241  and  9  is  19.89  degrees.  The  percentage  errors  are  2.3,  2.6  and  0.3  percent 
respectively. 

Another  example  in  Figure  9  and  Figure  8  shows  a  synthetic  dataset  and  its 
inversion  results.  The  model  is  a  spherical  object  located  in  a  constant  velocity 
medium.  The  synthetic  data  is  generated  in  3D  with  an  offset  of  approximately  five 
times  the  diameter  of  the  object.  The  imaging  is  done  with  GOCAD.  In  the  inversion 
image,  an  isosurface  is  created  for  the  peak  amplitude  so  that  the  spherical  object 
can  be  shown  in  3D  without  the  diffraction  events.  The  reflection  coefficient  at  the 
normal  incident  position  (the  top  of  the  sphere)  has  an  error  of  only  3.1  percent. 
The  position  around  the  top  of  the  sphere  is  comparably  accurate  except  in  the 
neighborhood  of  the  equator,  where  we  have  no  specular  returns.  There,  the  sphere  is 
filled  in  by  diffraction  returns  rather  than  by  true  specular  points  on  the  reflector.  The 
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Midpoint  (km) 

4  6  8  10 


Fig.  6.  Reflection  coefficient  R. 

lower  half  of  the  sphere  is  inaccurately  positioned  because  the  change  in  propagation 
speed  across  the  upper  sphere  surface  does  not  satisfy  the  criterion  of  being  a  depth 
dependent  background  velocity  model.  Thus,  1  used  a  constant  velocity  (1.5  km/s) 
down  to  the  water/seabed  interface  and  another  constant  below  (2.0  km/s).  Both  are 
lower  than  the  true  velocity  of  6  km/s  in  the  sphere.  Since  the  data  were  generated 
with  the  correct  velocity,  the  specular  returns  from  the  lower  part  of  the  sphere  came 
sooner  than  they  would  have  at  2.0  km/s.  Thus,  in  the  imaging,  the  lower  part  of 
the  sphere  is  pulled  up,  making  a  more  elliptical  surface.” 

CONCLUSIONS 

I  have  described  a  method  to  calculate  the  amplitude  weights  for  the  3-D  common 
offset  v{z)  inversion  by  ray  tracing.  The  traveltime  and  amplitude  tables  generated 
by  this  method  are  used  in  (1)  to  perform  Kirchhoff  inversion.  Tests  on  synthetic 
data  are  in  progress  to  improve  the  program  and  the  algorithm. 
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Fig.  7.  Synthetic  data. 
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Fig.  8.  Inversion  from  the  synthetic  data. 
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